home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / romberg.i < prev    next >
Text File  |  1995-07-26  |  4KB  |  121 lines
  1. /*
  2.    ROMBERG.I
  3.    Romberg integrator, after qromb in Numerical Recipes (Press, et.al.)
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func romberg(function, a, b, epsilon, notvector=)
  11. /* DOCUMENT integral= romberg(function, a, b)
  12.          or integral= romberg(function, a, b, epsilon)
  13.      returns the integral of FUNCTION(x) from A to B.  If EPSILON is
  14.      given, Simpson's rule is refined until that fractional accuracy
  15.      is obtained.  EPSILON defaults to 1.e-6.
  16.  
  17.      If the notvector= keyword is supplied and non-zero, then FUNCTION
  18.      may not be called with a list of x values to return a list of
  19.      results.  By default, FUNCTION is assumed to be a vector function.
  20.  
  21.      If the function is not very smooth, simpson may work better.
  22.  
  23.    SEE ALSO: simpson, max_doublings
  24.  */
  25. {
  26.   if (!epsilon || epsilon<0.) epsilon= 1.e-6;
  27.   a= double(a);
  28.   b= double(b);
  29.   s= array(0.0, 5);
  30.   h= [1.0, 0.25, 0.0625, 0.015625, 0.00390625];
  31.   for (i=1 ; i<=max_doublings ; ++i) {
  32.     ss= trapezoid(function, a, b, i, notvector);
  33.     if (i >= 5) {
  34.       s(5)= ss;
  35.       c= d= s;  /* Neville algorithm tableau */
  36.       ns= 4;    /* one less than smallest h, always last */
  37.       for (m=1 ; m<5 ; ++m) {
  38.     m5= 5-m;
  39.     ho= h(1:m5);
  40.     hp= h(m+1:5);
  41.     w= (c(2:m5+1)-d(1:m5))/(ho-hp);
  42.     d(1:m5)= hp*w;
  43.     c(1:m5)= ho*w;
  44.     dss= (2*ns<m5)? c(ns+1) : d(ns--);
  45.     ss+= dss;
  46.       }
  47.       if (abs(dss) < epsilon*abs(ss)) return ss;
  48.       /* extrapolation to h=0 always uses last 5 points */
  49.       s(1:4)= s(2:5);
  50.       h*= 0.25;
  51.     } else {
  52.       s(i)= ss;
  53.     }
  54.   }
  55.   error, "no convergence after "+pr1(2^i)+" function evaluations";
  56. }
  57.  
  58. func simpson(function, a, b, epsilon, notvector=)
  59. /* DOCUMENT integral= simpson(function, a, b)
  60.          or integral= simpson(function, a, b, epsilon)
  61.      returns the integral of FUNCTION(x) from A to B.  If EPSILON is
  62.      given, Simpson's rule is refined until that fractional accuracy
  63.      is obtained.  EPSILON defaults to 1.e-6.
  64.  
  65.      If the notvector= keyword is supplied and non-zero, then FUNCTION
  66.      may not be called with a list of x values to return a list of
  67.      results.  By default, FUNCTION is assumed to be a vector function.
  68.  
  69.      If the function is very smooth, romberg may work better.
  70.  
  71.    SEE ALSO: romberg, max_doublings
  72.  */
  73. {
  74.   if (!epsilon || epsilon<0.) epsilon= 1.e-6;
  75.   a= double(a);
  76.   b= double(b);
  77.   ost= os= -1.e-30;
  78.   for (i=1 ; i<=max_doublings ; ++i) {
  79.     st= trapezoid(function, a, b, i, notvector);
  80.     s= (4.*st - ost)/3.;
  81.     if (abs(s-os) < epsilon*abs(os)) return s;
  82.     os= s;
  83.     ost= st;
  84.   }
  85.   error, "no convergence after "+pr1(2^i)+" function evaluations";
  86. }
  87.  
  88. local max_doublings;
  89. /* DOCUMENT max_doublings= 20
  90.      is the maximum number of times romberg or simpson will split the
  91.      integration interval in half.  Default is 20.
  92.  */
  93. max_doublings= 20;
  94.  
  95. func trapezoid(function, a, b, n, notvector)
  96. {
  97.   extern _trapezoid_i, _trapezoid_s;
  98.   if (n==1) {
  99.     _trapezoid_i= 1;
  100.     return _trapezoid_s= 0.5*(b-a)*(function(a)+function(b));
  101.   }
  102.   dx= (b-a)/_trapezoid_i;
  103.   if (!notvector) {
  104.     /* conserve memory-- dont try more than 8192 points at a time */
  105.     if (_trapezoid_i <= 8192) {
  106.       s= sum(function(span(a,b,_trapezoid_i+1)(zcen)));
  107.     } else {
  108.       x= a+(indgen(8192)-0.5)*dx;
  109.       s= sum(function(x));
  110.       dx2= 8192*dx;
  111.       for (i=8193 ; i<=_trapezoid_i ; i+=8192) s+= sum(function((x+=dx2)));
  112.     }
  113.   } else {
  114.     x= a+0.5*dx;
  115.     s= function(x);
  116.     for (i=2 ; i<=_trapezoid_i ; ++i) s+= function((x+=dx));
  117.   }
  118.   _trapezoid_i*= 2;
  119.   return _trapezoid_s= 0.5*(_trapezoid_s + dx*s);
  120. }
  121.